# load data
dat <- read_xlsx(here("Data", "Hodkinson_2020.xlsx"))

# preprocess
dat$vi <- dat$`Std Err2`^2
dat <- dat[, c("Intervention model", "Outcome", "Hedgesg", "vi", "Study name")]
names(dat) <- c("dimension1", "dimension2", "yi", "vi", "study_id")

# show data in a datatable
t1 <- dat %>%
    dfround(3) %>%
    DT::datatable()
t1

# get estimate for each cell
est_dat <- dat %>% group_by(dimension1, dimension2) %>%
  group_modify(~ meta_aggregate(.x, rho = 0.5)) %>% ungroup()

# traditional map with the number of study
est_dat$dimension1 <- as.factor(est_dat$dimension1)

Box1_map1 <- ggplot(est_dat, aes(x = dimension1, y = dimension2, size = n_studies) ) +
  geom_point(alpha = 0.5, color = RColorBrewer::brewer.pal(n = 8, name = "Dark2")[1]) + 
  labs(x = "Dimension 1 (Intervention)", y = "Dimension 2 (Outcome)") +
  scale_size(range=c(5,10)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) + # 
  scale_y_discrete(labels = function(x) str_wrap(x, width = 10)) +
  theme_bw() +
  guides(size = "none") +
  theme(legend.position='top', 
        legend.justification='right',
        legend.direction='horizontal', 
        axis.text = element_text(color = "black"),
        axis.title = element_text(color = "black")) +
  geom_text(aes(label = as.character(n_studies)), size = 4, color = "gray10") +
  labs(caption = "The value in the cell is the number of studies") + 
   theme(plot.caption = element_text(size = 10, color = "gray10", face = "italic"),
         axis.text = element_text(size = 12),
         axis.title = element_text(size = 12, face = "bold")) 


# traditional map with the number of effect size
Box1_map2 <- ggplot(est_dat, aes(x = dimension1, y = dimension2, size = n_es)) +
  geom_point(alpha = 0.5, color = RColorBrewer::brewer.pal(n = 8, name = "Dark2")[3]) + 
  labs(x = "Dimension 1 (Intervention)", y = "Dimension 2 (Outcome)") +
  scale_size(range=c(5,10)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  scale_y_discrete(labels = function(x) str_wrap(x, width = 10)) +
  theme_bw() +
  guides(size = "none") +
  theme(legend.position='top', 
        legend.justification='right',
        legend.direction='horizontal',  
        axis.text = element_text(color = "black"),
        axis.title = element_text(color = "black")) +
  geom_text(aes(label = as.character(n_es)), size = 4, color = "gray10") +
  labs(caption = "The value in the cell is the number of effect sizes") +
   theme(plot.caption = element_text(size = 10, color = "gray10", face = "italic"),
         axis.text = element_text(size = 12),
         axis.title = element_text(face = "bold", size = 12))



# with effect size information
Box1_map3 <- ggplot(est_dat, aes(x = dimension1, y = dimension2, size = n_es, color = estimate)) +
  geom_point(alpha = 0.6) + 
  scale_color_gradient(
    low = "blue",
    high = "red", 
    limits = c(-1,1),
    guide = "colourbar") + 
  labs(x = "Dimension 1 (Intervention)", y = "Dimension 2 (Outcome)", color = "Meta-analytic mean effect size") +
  scale_size(range=c(5,10)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  scale_y_discrete(labels = function(x) str_wrap(x, width = 10)) +
  theme_bw() +
  guides(size = "none") +
  theme(legend.position='top', 
        legend.justification='right',
        legend.direction='horizontal', 
        axis.text = element_text(color = "black"),
        axis.title = element_text(color = "black")) +
  geom_text(aes(label = as.character(n_es)), size = 4, color = "gray10") +
  labs(caption = "The value in the cell is the number of effect sizes") + 
   theme(plot.caption = element_text(size = 10, color = "gray10", face = "italic"),
         axis.text = element_text(size = 12),
         axis.title = element_text(face = "bold", size = 12))


# 3-dimensional graph - we do not recommend to 3 dimensional graph because this is not that friendly to end-user
dat <- read_xlsx(here("Data","Hodkinson_2020.xlsx"))
dat$vi <- dat$`Std Err2`^2
dat <- dat[, c("Intervention model", "Outcome", "Age group", "Hedgesg", "vi", "Study name")]
names(dat) <- c("dimension1", "dimension2", "dimension3", "yi", "vi", "study_id")

# get estimate each cell 
est_dat <- dat %>% group_by(dimension1, dimension2, dimension3) %>%
  group_modify(~ meta_aggregate(.x, rho = 0.5)) %>% ungroup()

# with the third dimension
Box1_map4 <- ggplot(est_dat, aes(x = dimension1, y = dimension2, 
                      shape = dimension3,
                      size = n_es  
                      )) +
  geom_point(alpha = 0.6,  color = RColorBrewer::brewer.pal(n = 8, name = "Dark2")[5], aes(group = dimension3), position = position_dodge(width= 0.5)) +
 scale_color_gradient(
    low = "blue",
    high = "red", 
    limits = c(-1,1),
    guide = "colourbar") + 
  scale_shape_manual(values = c(15, 19, 17, 18, 10, 4, 3, 7, 8, 13)) +
  labs(x = "Dimension 1 (Intervention)", y = "Dimension 2 (Outcome)", shape = "Diamension 3 (Population)") +
  scale_size(range=c(5,10)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  scale_y_discrete(labels = function(x) str_wrap(x, width = 10)) +
  theme_bw() +
  guides(size = "none") +
  theme(legend.position='top', 
        legend.justification='right',
        legend.direction='horizontal', 
        axis.text = element_text(color = "black"),
        axis.title = element_text(color = "black")) +
  geom_text(aes(label = as.character(n_es), group = dimension3), size = 4, color = "gray10", position = position_dodge(width = .5)) +
  labs(caption = "The value in the cell is the number of effect sizes") + 
   theme(plot.caption = element_text(size = 10, color = "gray10", face = "italic"),
         axis.text = element_text(size = 12),
         axis.title = element_text(face = "bold", size = 12))


#Box1_map1 + Box1_map2 + Box1_map3 + Box1_map4 + patchwork::plot_layout(ncol = 2, nrow = 2) + patchwork::plot_annotation(tag_levels = "A")
#png(filename = "Box1_map.png", width = 10, height = 10, units = "in", type = "windows", res = 400)
plot_grid(Box1_map1, Box1_map2, Box1_map3, Box1_map4, labels = c('A','B','C','D'), label_size = 14, nrow = 2, ncol = 2)

# load data
dat <- read.csv(here("Data", "Mertens_2021.csv"))
names(dat)[1] <- "publication_id"
# show data in a data
t2 <- dat %>%
    dfround(3) %>%
    DT::datatable()
t2

# pre-process
dat <- dat %>%
    dplyr::select(type_experiment, intervention_technique, intervention_category,
        domain, population, location)

# png(filename = './Box1_sankey.png', width = 6, height = 5, units = 'in', type
# = 'windows', res = 400)
ggplot(dlong(dat, type_experiment, location, population, intervention_category, domain,
    intervention_technique), aes(x = x, next_x = next_x, node = node, next_node = next_node,
    fill = factor(node), label = node)) + sankey_p(flow.alpha = 0.8, node.color = "transparent") +
    sankey_p_label(size = 3, color = "white", fill = "gray10", alpha = 0.6) + ggsci::scale_color_tron() +
    theme_sankey(base_size = 10) + labs(x = NULL) + theme(legend.position = "none",
    plot.title = element_text(hjust = 0.5), axis.text.x = element_text(color = "black",
        size = 9)) + scale_x_discrete(labels = c("Moderator 1 \n (Experimental \n approach)",
    "Moderator 2 \n (Geographical \n location)", "Moderator 3 \n (Population \n characteristic)",
    "Moderator 4 \n (Architecture \n category)", "Moderator 5 \n (Behavioral \n domain)",
    "Moderator 6 \n (Intervention \n technique)"), position = "top")

# load data
dat <- read_xlsx(here("Data", "Sanders_2021.xlsx"))
# show data in a data
t3 <- dat %>%
    dfround(3) %>%
    DT::datatable()
t3

# some data wrangling
dat$Species <- gsub("\\.", "_", dat$Species)
dat <- dat[!is.na(dat$Species), ]  # 

# check species we have length(unique(dat$Species)) #184 unique species names,
# if no misspelling

# find Open Tree Taxonomy (OTT) IDs for each species
taxa <- tnrs_match_names(names = unique(dat$Species))

# for illustrative purpose, we just delete the species that are not matched
# find location
pos <- which(dat$Species %in% c("Phodopus_sungeroru", "Electra?__sp", "Molgula_sp",
    "Baetis_spp", "Agrostis_tenuis", "Anthoxanthum_odoratum", "Myotis_daubentonii",
    "_Myotis_mystacinus", "_Myotis_brandtii", "_Myotis_nattereri", "_Plecotus_auritus",
    "Pipistrellus_pipistrelles", "_Pipistrellus_nathusii", "Nyctalus_nyctalus", "_Nyctalus_leisleri",
    "_Eptesicus_serotinus_", "All_families", "Pipistrellus_hesperidus/Hypsugo_anchietaia",
    "Myotis?_sp", "MCF7", "7288cct", "A__cristatellus_", "A__evermanni", "A__gundlachi",
    "A__sagrei", "bat_community_", "N__leisleri", "Poecile_sp", "Coleoptera", "_Diptera",
    "_Lepidoptera", "_Erebidae", "_Chironomidae", "_Noctuidae_and_Psychodidae", "Nyctalus_and_Eptesicus_spp",
    "Anoectochilus_roxburghii", "Myotis_pilosatibialis", "Silene_latifolia"))
dat <- dat[-pos, ]

# match again
taxa <- tnrs_match_names(names = unique(dat$Species))  # still some are not matched. But let's resolve it latter

# check whether occur in the synthetic tree
ott_in_tree <- ott_id(taxa)[is_in_tree(ott_id(taxa))]
# length(ott_id(taxa)) - length(is.na(ott_in_tree)) # still 10 did not appear
# in synthetic tree

# check which 10
out_tree <- filter(taxa, !ott_id %in% ott_in_tree)
# out_tree$search_string

dat <- dat[!dat$Species %in% c("Leptocythere_pellucida", "Myotis_daubentonii,_Myotis_mystacinus,_Myotis_brandtii,_Myotis_nattereri,_Plecotus_auritus",
    "Pipistrellus_pipistrelles,_Pipistrellus_nathusii", "Nyctalus_nyctalus,_Nyctalus_leisleri,_Eptesicus_serotinus_",
    "Myotis_lucifugus", "Acyrthosiphont", "Carabidae", "Staphylinidae", "Eptesicus_bottae",
    "Corethrella_spp", "P__pipistrellus", "Myotis", "Coleoptera,_Diptera,_Lepidoptera,_Erebidae,_Chironomidae,_Noctuidae_and_Psychodidae",
    "Pipistrellus"), ]

# match again
taxa <- tnrs_match_names(names = unique(dat$Species))
# check again whether all otts occur in the synthetic tree
ott_in_tree <- ott_id(taxa)[is_in_tree(ott_id(taxa))]
# length(ott_id(taxa)) - length(is.na(ott_in_tree)) # all good

# now every ott occur in the synthesistic tree. But for the sake of brevity, we
# only visualize a subset of the tree only use studies measuring physiology
Physiology <- subset(dat, Category == "Physiology")
# remove NA
Physiology <- subset(Physiology, Species != "NA")
taxa <- tnrs_match_names(names = unique(Physiology$Species))

# make phylo tree
tree <- suppressWarnings(tol_induced_subtree(ott_ids = ott_id(taxa)))

# the tip labels contain OTTs, which means they will not perfectly match the
# species names in our dataset or the taxon map that we created earlier, remove
# the extra information from the tip labels later; with the IDs removed, we can
# use our taxon map to replace the tip labels in the tree with the species
# names from dataset
tree$tip.label <- strip_ott_ids(tree$tip.label, remove_underscores = TRUE)

# decapitalise species names to match with the search string names in taxa
Physiology <- Physiology %>%
    mutate(search_string = tolower(Species))

# align data
Physiology <- left_join(Physiology, dplyr::select(taxa, search_string, unique_name,
    ott_id), by = "search_string")

# create the variables of spp and phylo
Physiology <- Physiology %>%
    mutate(spp = unique_name, phylo = unique_name)

# prepare annotation data species-specific estimates as fixed effect
agg.es <- escalc(measure = "SMD", m1i = Experimental_Mean, m2i = Control_Mean, sd1i = Experimental_SD,
    sd2i = Control_SD, n1i = Experimental_N, n2i = Control_N, data = Physiology) %>%
    aggregate(cluster = phylo, struct = "CS", rho = 0.5, addk = TRUE)
agg.es <- agg.es[c("phylo", "yi", "vi")]
names(agg.es) <- c("Species", "Mean", "SE")
agg.es <- agg.es %>%
    mutate(Lower_bound = Mean - sqrt(SE) * qnorm(0.975), Upper_bound = Mean + sqrt(SE) *
        qnorm(0.975)) %>%
    arrange(Species)
N_obs <- Physiology %>%
    group_by(phylo) %>%
    summarise(N_obs = n())  # calculate sample size
names(N_obs)[1] <- "Species"
fe.spp.es <- left_join(agg.es, N_obs, by = "Species")
tip.label <- data.frame(Species = tree$tip.label)  # extract tip label
fe.spp.es2 <- left_join(tip.label, fe.spp.es, by = "Species")
fe.spp.es2 <- fe.spp.es2 %>%
    mutate(z = Mean/SE, p = pnorm(abs(z), lower.tail = F) * 2)

# get class data
spp.class <- dplyr::distinct(Physiology, phylo, .keep_all = TRUE) %>%
    dplyr::select(Class, phylo)
names(spp.class)[2] <- "Species"
spp.class2 <- left_join(tip.label, spp.class, by = "Species")

# make each panel
tree.p1 <- ggtree(tree, layout = "rectangular", cex = 0.4)
tree.p2 <- tree.p1 %<+% spp.class2 + geom_tiplab(aes(color = Class), size = 3, fontface = "italic",
    align = T, offset = 0.3) + geom_tippoint(aes(color = Class)) + xlim_expand(c(0,
    21), panel = "Tree") + guides(color = "none") + scale_color_viridis_d()
tree.p3 <- tree.p2 + geom_facet(panel = "Effect size", data = fe.spp.es2, geom = ggstance::geom_pointrangeh,
    mapping = aes(x = Mean, xmin = Lower_bound, xmax = Upper_bound, color = Class)) +
    theme_tree2() + theme(strip.background = element_rect(fill = "white"))

# png(filename = 'Box1_phylo.png', width = 6, height = 5, units = 'in', type =
# 'windows', res = 400)
facet_widths(tree.p3, c(Tree = 0.4, `Effect size` = 0.2))

# load data
dat <- read_xlsx(here("Data", "Sanders_2021.xlsx"))
# load bibliographic data
M <- convert2df(here("Data", "Sanders_2021_bib.csv"), dbsource = "scopus", format = "csv") %>%
    suppressMessages()  # note that using here function does work for data wrangling

# merge data
dat.M <- left_join(M, dat, by = c(DI = "DOI"))  # DI denote DOI in bib data

# construct author collaboration network
NetMatrix <- biblioNetwork(M, analysis = "collaboration",  network = "authors", sep = ";")
net_matrix <- as.matrix(NetMatrix)
diag(net_matrix) <- 0 
g <- graph_from_adjacency_matrix(net_matrix, mode = "lower", weighted = "weight")
# computing centrality measures for each vertex 
V(g)$indegree <-  igraph::degree(g, mode = "in")
V(g)$outdegree <- igraph::degree(g, mode = "out")
V(g)$closeness <- igraph::closeness(g, mode = "total")
V(g)$betweeness <- igraph::betweenness(g, normalized = TRUE)

#  graph
set.seed(2023)
# using random walks to detect the network community
wtc <- cluster_walktrap(g)
# modularity(wtc)
# modularity(g, membership(wtc))
# member <- communities(wtc)
# n_cluster <- sapply(member, length)

#png(filename = "./Box2_author.png", width = 15, height = 15, units = "in", type = "windows", res = 400)
plot(wtc, g, 
     vertex.size = V(g)$indegree/3, # or V(g)$outdegree + 1
     vertex.label = NA,
     edge.arrow.size = .25)

# vertex and edges attributes can be accessed via the V and E functions,
# respectively we can list what vertex/edge attributes are available:
# list.vertex.attributes(g) list.edge.attributes(g) V(g)$name[1:6]

# extracting each vectex features as a data.frame
stats <- as_data_frame(g, what = "vertices")

# computing quantiles for each vertex attributes
stats_degree <- with(stats, {
    cbind(indegree = quantile(indegree, c(0.025, 0.5, 0.975), na.rm = TRUE), outdegree = quantile(outdegree,
        c(0.025, 0.5, 0.975), na.rm = TRUE), closeness = quantile(closeness, c(0.025,
        0.5, 0.975), na.rm = TRUE), betweeness = quantile(betweeness, c(0.025, 0.5,
        0.975), na.rm = TRUE))
})

# compute some statistics at the graph level:
cbind(size = vcount(g), nedges = ecount(g), density = edge_density(g), recip = reciprocity(g),
    centr = centr_betw(g)$centralization, pathLen = mean_distance(g))

# country
M_country <- metaTagExtraction(M, Field = "AU_CO", sep = ";")
NetMatrix <- biblioNetwork(M_country, analysis = "coupling", network = "countries",
    sep = ";")

net_matrix <- as.matrix(NetMatrix)
# get rid of collaboration with same country
diag(net_matrix) <- 0
# getting rid of lower triangle
net_matrix[lower.tri(net_matrix)] <- 0

# colnames(net_matrix) - change to title case:
colnames(net_matrix) <- str_to_title(colnames(net_matrix))
# rownames(net_matrix) - change to title case:
rownames(net_matrix) <- str_to_title(rownames(net_matrix))
# Fix 'Usa' to 'United States' :
colnames(net_matrix)[colnames(net_matrix) == "Usa"] <- "United States"
rownames(net_matrix)[rownames(net_matrix) == "Usa"] <- "United States"

# color palette
color <- viridis::viridis(34, alpha = 1, option = "D")
color <- color[sample(1:34)] 

circos.clear()
circos.par(start.degree = 90, 
           gap.degree = 3,  
           points.overflow.warning = FALSE)

#png(filename = "Box2_country.png", width = 5, height = 5, units = "in", type = "windows", res = 400)
chordDiagram(net_matrix, 
             grid.col = color, transparency = 0.1, 
             directional = 1,direction.type = c("arrows", "diffHeight"), diffHeight  = -0.04, annotationTrack = "grid", annotationTrackHeight = c(0.05, 0.1),link.arr.type = "big.arrow",link.sort = TRUE, link.largest.ontop = TRUE,
             preAllocateTracks = 1)
# add text and axis
circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
  xlim = get.cell.meta.data("xlim")
  ylim = get.cell.meta.data("ylim")
  sector.index = get.cell.meta.data("sector.index")
  # add names to the sector. 
  circos.text(x = mean(xlim), y = ylim[1] + .01, 
              labels = sector.index, facing = "clockwise", # clockwise
              niceFacing = TRUE, cex = 0.6, adj = c(0, 0.5))
  # add graduation on axis
  circos.axis(h = "top", labels.cex = 0.01, major.tick.length = 0.01, sector.index = sector.index, track.index = 2, labels.niceFacing = T)
}, bg.border = NA)